%% The actual Travel model and modifications
%
% Input:  PreferenceFactor (used for minimization),cfg,exvars,depvars,indepvars,r,t
% Output: Residuals (depeding on calibration type; only during calibrations)
%         depvars structure
%
function varargout = ModelGcam(cfg,exvars,depvars,indepvars,r,time,varargin)
%% Destinguish between calibration or validation and specific type of run
for n = 2:2:length(varargin)
    eval( [varargin{n-1} '(r,1,:)=varargin{n};']);
end
modes  =  indepvars.UseModes(r,1,:); % get rid of modes not used by Schafer (e.g. High speed train)

exvars.Speed = repmat(exvars.Speed(:,35,:),[1,130,1]);
if isfield(varargin,'indepvars.Sigma') &&1==0
        indepvars.Sigma(r)  = ...
        sum(exvars.TravelDemand(r,cfg.CalibrationYears(end),modes),3)...
        ./((exvars.Income(r,cfg.CalibrationYears(end))^indepvars.IncomeElasticity(r)).* ...
        (sum((exvars.TravelPrice (r,cfg.CalibrationYears(end),modes)...
        + exvars.WageRate (r,cfg.CalibrationYears(end))...
        ./ exvars.Speed(r,cfg.CalibrationYears(end),modes))...
        .*exvars.ModeSplit(r,cfg.CalibrationYears(end),modes),3)^-indepvars.PriceElasticity(r)));
end


if cfg.ConvergencePFandSW && sum(r==[1 2 3])==0
    depvars.ShareWeight(r,:,:) = repmat(indepvars.ShareWeight(r,1,:),[1 size(depvars.TravelDemand,2) 1]);
    depvars.ShareWeight(r,:,:) = repmat(indepvars.ShareWeight(r,1,:),[1 size(depvars.TravelDemand,2) 1]);
    depvars.ShareWeight(r,:,3) = depvars.ShareWeight(r,:,4).*((exvars.Income(1,:).^0.5)./(exvars.Income(r,:).^0.5)).*(indepvars.SWSaturation(r).^-0.5);    
else
    depvars.ShareWeight(r,:,:) = repmat(indepvars.ShareWeight(r,1,:),[1 size(depvars.TravelDemand,2) 1]);
end

for t=time
%% Gcam model
    depvars.ModeCostGcam(r,t,modes)      ...
        = indepvars.TravelPrice (r,t,modes) + exvars.WageRate (r,t)./ exvars.Speed(r,t,modes);

    
    depvars.ModeSplit(r,t,modes)         ...
        = ((depvars.ShareWeight(r,t,modes) .* depvars.ModeCostGcam (r,t,modes)).^-indepvars.CostDistrParameter(r))...
        ./ (sum((depvars.ShareWeight(r,t,modes) .* depvars.ModeCostGcam(r,t,modes)).^-indepvars.CostDistrParameter(r),3));

    depvars.TravelDemand (r,t,modes)     ...
        = (indepvars.Sigma(r)*(exvars.Income(r,t)^indepvars.IncomeElasticity(r)).* ...
        (sum(depvars.ModeCostGcam(r,t,modes).*depvars.ModeSplit(r,t,modes),3)^-indepvars.PriceElasticity(r))).* ...
        depvars.ModeSplit(r,t,modes);
%% Errors
    if t<36 % ensures function is only called with historic data
        depvars.ModeWeight(r,t,modes) = sqrt(exvars.TravelDemand(r,t,modes))./sum(sqrt(exvars.TravelDemand(r,t,modes)),3);% Ensures not all modes are equal, but that car is not the only important mode
        depvars.ObjectiveFunction(r,t) = sqrt(sum((depvars.TravelDemand(r,t,modes) - exvars.TravelDemand(r,t,modes)).^2 ...
        ./ exvars.TravelDemand(r,t,modes).^2 .* depvars.ModeWeight(r,t,modes).^2,3)); % used for optimization
        depvars.Resid(r,t)             =  depvars.ObjectiveFunction(r,t);
        
        depvars.errors.TDfit(r,1)    = sqrt(sum(((sum(depvars.TravelDemand(r,cfg.CalibrationTime,modes),3) - sum(exvars.TravelDemand(r,cfg.CalibrationTime,modes),3)) ...
    ./sum(exvars.TravelDemand(r,cfg.CalibrationTime,modes),3)).^2,2)/length(cfg.CalibrationTime));

     depvars.errors.ModeFit(r,1)  = sqrt(sum(sum((exvars.TravelDemand(r,cfg.CalibrationTime,modes)...
         ./repmat(sum(exvars.TravelDemand(r,cfg.CalibrationTime,modes),3),[1 1 length(cfg.modes(modes))])-depvars.ModeSplit(r,cfg.CalibrationTime,modes)).^2,3))/length(cfg.CalibrationTime)); 
    end
end

%% Write output
if isfield(varargin,'indepvars.Sigma')
    varargout{1} = mean(depvars.errors.TDfit(r,:));
elseif isfield(varargin,'indepvars.ShareWeight')&& ~cfg.OFout
    varargout{1} = mean(depvars.errors.ModeFit(r,:));
else
    varargout{1} = mean(depvars.ObjectiveFunction(r,time));
end
    
if nargout>=2
    varargout{2} = depvars;
end
if nargout >=3
    varargout{3} = indepvars;
end

%% Visualizes minimization process
if cfg.LiveMinimizeVisuals{1} > 1 && ismember(time(1),cfg.LiveMinimizeVisuals{2}); % bit clumsy to make the code faster
    GeneralPlotting('LiveMinimizeVisuals','Dont save',cfg,exvars,depvars,indepvars,r,time);
elseif cfg.PlotMinimize == 1
    GeneralPlotting('SaveToGlobal',' ',cfg,exvars,depvars,indepvars,r,time);
end
    
end